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Abstract: We address the issue of computing the non-linear matter power spec- 
trum on mildly non-linear scales with efficient semi-analytic methods. We implemented 
M. Pietroni's Time Renormalization Group (TRG) method and its Dynamical 1-Loop 
(D1L) limit in a numerical module for the new Boltzmann code CLASS. Our publicly released 
module is valid for ACDM models, and optimized in such a way to run in less than a minute 
for D1L, or in one hour (divided by number of nodes) for TRG. A careful comparison of the 
D1L, TRG and Standard 1-Loop approaches reveals that results depend crucially on the 
assumed initial bispectrum at high redshift. When starting from a common assumption, 
the three methods give roughly the same results, showing that the partial resumation of 
diagrams beyond one loop in the TRG method improves one-loop results by a negligible 
amount. A comparison with highly accurate simulations by M. Sato & T. Matsubara shows 
that all three methods tend to over-predict non-linear corrections by the same amount on 
small wavelengths. Percent precision is achieved until k ~ 0.2/iMpc -1 for z > 2, or until 
k ~ 0.14 /iMpc" 1 at z = 1. 
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1. Motivations 

Large scale structures in our universe have formed during the matter dominated era, start- 
ing from a very homogeneous state. In order to explain this mechanism, one must compute 
the evolution of the dominating species during this period: Cold Dark Matter (CDM) and 
baryons in the standard ACDM Model. It is generally agreed that structures formed from 
the gravitational collapse of small density perturbations. Depending on their wavelength, 
they entered at different times within the Hubble scale, which plays the role of a causal 
horizon for this process. Perturbations with very large wavelength, or very small wave 
number, had no time to evolve significantly beyond the linear regime. In practice, one can 
apply the linear perturbation theory for wave numbers smaller than 0.1/iMpc -1 . In the op- 
posite limit, for scales under ~ 10 Mpc today, the presence of highly non-linear structures 
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such as galaxies points out the need for completely non-linear computations. This is taken 
care of by N-body simulations, that follow the evolution of a large ensemble of galaxies. 

Current and upcoming surveys such as the Sloan Digital Sky Survey , the Large Syn- 
optic Survey Telescope 2 and other Large Scale Structure (LSS) experiments probe this 
evolution with increasingly high precision. Through these new observations, cosmology 
opens a new window to confirm, constrain or infirm different high energy physics scenarios, 
related for instance to: neutrino masses, inflationary and dark energy models, or modifica- 
tions of gravity. 

The discriminating power of LSS observations depends crucially on the maximum 
wavenumber used in the comparison with the theory. By limiting the analysis to linear 
scales with k < k max = 0.1/iMpc -1 , one loses a lot of sensitivity, since the total amount of 
information scales like fc^az- For instance, the strong dependence of neutrino mass error 
bars on k max is illustrated in Q . To further enhance the sensitivity, one would be tempted 
to use highly detailed N-body simulations. Unfortunately, in order to find both the best- 
fitting values and the error bars of the free cosmological parameters of a given scenario, one 
needs to compute a huge number of theoretical spectra corresponding to different points 
in parameter space. With the most efficient techniques (Monte Carlo Markov Chains), 
a minimum of lO'OOO to lOO'OOO points is necessary, depending on the complexity of the 
model. N-body simulations are far too slow for being carried in each of these points, or 
even in sizable fraction of them. This raises the need for semi-analytical tools to make 
accurate predictions on interesting scales. Even if such tools remain accurate only in a 
small range of mildly non-linear scales (just above 0.1/iMpc -1 ), they can play a crucial role 
in measuring quantities like neutrino masses. 

Some of these semi-analytical tools just consist in fitting formulas calibrated to N- 
body simulation, as for instance HALOFIT Q. HALOFIT represents the fastest thinkable 
way to account for non-linear perturbations, but its range of validity is limited to the 
minimal ACDM model. Extending it to non-minimal cosmological scenarios requires many 
N-body simulations (and in some cases, like for models containing species with a sizable 
free-streaming length, carrying any N-body simulation remains very involved). Moreover, 
HALOFIT does not provide a good fit to the Baryon Acoustic Oscillations (BAOs) in the 
non-linear spectrum. 

Other semi-analytical methods have been proposed to actually calculate the non-linear 
power spectrum in Fourier space, taking into account the effects of mode coupling to some 
extent (for reviews and comparisons, see e.g. ||, ||, ||, ||, |7j §). Of course, all these 
approaches fail when dealing with the highly non-linear regime. However, their formulation 
stays consistent within the mildly non-linear regime, up to some k max which depends on 
the method. Any tool implementing a good compromise between computing time on the 
one hand, and accuracy (i.e. large k max ) on the other hand, can be extremely useful for two 
purposes: first, for calibrating fitting formulas in the mildly non- linear regime for extended 
cosmological scenarios, without running N-body simulations; second, if the method is fast 

x http : //www. sdss . org 
2 http : //www. lsst . org 
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enough, for being employed at each point in parameter space when fitting cosmological 
models to the data. 



In this paper, we concentrate on the Time Renormalization Group (TRG) method 
proposed by Massimo Pietroni . This method consists in integrating over time a coupled 
system of differential equations for the density and velocity power spectra and bispectra. 
Since the non-linearity computed at a given time step affect the time-derivative of the 
spectra at this step, the TRG method continuously includes higher-order corrections to the 
linear power spectra, and can be seen as a simple way to resum a sub-class of diagrams 
beyond one loop. 

A simple variant of the TRG equations, consisting in using products of the linear 
power spectra in the non-linear source term of the equations, provides strictly one-loop 
results. However, the one-loop power spectrum computed in that way is obtained from 
dynamical equations (we actually call this method D1L for Dynamical 1-Loop). This 
should be contrasted with the Standard 1-Loop (S1L) method in which the time evolution 
is integrated away, using some simplifying assumptions. The main focus of this paper 
consists in a detailed comparison between S1L, D1L and TRG results for realistic ACDM 
models. We will see that assumptions concerning the initial bispectrum at high redshift 
crucially affect the results al low redshift. This point had been overlooked in the past, and 
will lead us to the conclusion that when starting from the same assumptions, the three 
methods only differ by a negligible amount. Hence, using the TRG method (at the order 
discussed in the current literature) instead of the much faster D1L algorithm does not 
appear to be justified. 

We also provide some details about our numerical implementation of these methods 
in the form of a C module for the new Boltzmann code CLASS 3 |9], 10[, released with the 



version vl.2 of the code. Our implementation allows for a small computing time (for the 
TRG method, 70 minutes on a single-core CPU, to be divided by the number of cores in 
a parallel execution; and for D1L, less than a minute), opening the possibility to integrate 
it within a parameter extraction code. These approaches could easily be used e.g. for 
computing CMB lensing or cosmic shear power spectra in the mildly non-linear regime, 
without recurring to any N-body simulation or empirical fitting formulas. They also provide 
the non-linear velocity and cross density-velocity power spectra. 

We review the standard one-loop approach very briefly in section ^, M. Pietroni's TRG 
method in section |3|, and our numerical implementation of the previous equations in sec- 
tion ^. We present some self-consistency checks in the exact Einstein-De-Sitter (EdS) limit 
in [|, and show our results for a realistic ACDM model in section ||. We include a comparison 
with N-body simulations performed with GADGET [|ll]] by Sato & Matsubara ||, and with 
HALDFIT. A summary is provided in section [?]. Appendix |A| recalls how the continuity and 
Euler equations are derived from the Boltzmann equation. Finally, Appendix [B| provides 
details on the code itself. 



? http : / / class-code . net 



- 3 - 



2. Standard one-loop approach 



The one-loop density power spectrum is usually computed at any redshift using the formula 

PNL(k) = P L (k) (2.1) 
, 2tt f°° . , / f k+ P , ( {2k 4 + 3k 2 (q 2 +p 2 ) + lOpV - 5(p 4 + <z 4 )) 2 \ 

+ T y o ^l(p) ^ ^ j 

fc 3 P L (fc) / 12A; 2 100p 2 42p 4 3fc 3 (7p 2 M0 \fp 2 ,\\ I 1 ~ k\ \ \ 



+ 



where Ph{k) is the linear density power spectrum (see e.g. || and references therein). This 
result, that we will call S1L for Standard 1-Loop, is obtained by expanding the density and 
velocity fields describing a single-flow pressureless fluid at order two in perturbations. It 
is based however on additional approximations concerning the time-evolution of linear and 
non-linear perturbations. 

Indeed, the above formula corresponds to the exact one- loop solution of the continuity 
and Euler equations for such a fluid only in the perfect Einstein-De-Sitter limit (more pre- 
cisely, assuming that all modes obey to a newtonian evolution in a fully matter-dominated 
universe). In this case, the matter density perturbations can be expanded as 

oo 

6(r,k) = Y j a(r) i 5^(k) (2.2) 

i=i 

where r stands for conformal time, a % for the scale factor to the power i, and 5^' for the 
i-th order perturbation rescaled by a*. A similar expansion is performed for the velocity 
divergence field 9: 

oo 

e(r,k) = H(T)J2^rye^(k) , (2.3) 

with % = a! I a. When writing the relations between the 6® and 0® functions, all time- 
dependent factors simplify, thanks to the structure of the non-linear continuity and Euler 
equations. After collecting all terms, one can write the non- linear power spectrum as a 
function of the linear one evaluated at the same redshift, without any integral over time. 

In a realistic ACDM universe, the linear density fluctuations undergo a complicated 
evolution at high redshift: their behavior is radically different on super /sub-Hubble scales, 
during radiation/matter domination, before/after photon decoupling. After recombination, 
when radiation perturbations become negligible, their evolution can be caught by a scale- 
independent growth factor often parametrized as a(r) x g(r), where g is a function of time 
falling below one close to A or Dark Energy domination. Similarly, the linear velocity 
divergence 9 scales like a(r) x /(r), with / = g + ag'/a'. We show g(z) and f(z) for typical 
ACDM parameters in figure |]. As long as one can neglect massive neutrinos or eventual 
dark energy perturbations, high-order terms in perturbation theory should still be separable 
functions of time and wavenumber, because the model contains no characteristic length that 
could lead to a scale-dependent growth factor within the Hubble radius. However, these 
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growth factor vs. redshift 




Figure 1: Functions g(z) and f(z) related to the linear growth factor of density perturbations and 
velocity divergence, for a ACDM model with Qa = 0.735. 

terms do not scale like powers of [a g] . This can be checked by trying an expansion of the 
type 

oo oo 

5(t, k) = ^Hr)g(r)Y5^(k) , 9{r, k) = U(t) ^[a^g^fe^ik). (2.4) 

i=l i=l 

Since the non-linear equations of motion still contain factors in a'/a, not [ag]'/[ag], the 
time-dependence cannot be factored out in the same way as in the Einstein-De-Sitter case 
(for more details see e.g. [||, |6|). In summary, the formula fl2.1| ) is only approximate, due 
to the non-trivial evolution of perturbation in a realistic ACDM universe, both at high 
redshift (recombination time and before) and at low redshift (A domination). 

In order to evaluate the corresponding error, one should compare the power spectrum 
obtained from eq. fl2.1|) with that obtained from a full integration over time of one-loop 
equations. This task becomes easy once the TRG equations have been correctly imple- 
mented, since the TRG method relies on dynamical equations, and since the one-loop 
order corresponds to a simple limit of these equations. A test presented in § suggests that 
the error induced by low-redshift effects is very small. Below, we will confirm this finding, 
but interestingly, we will show that the error induced by mistreating the initial bispectrum 
at high redshift is instead very significant. 

3. Time Renormalization Group equations 

For the rest of this article, the following conventions are assumed. The universe is con- 
sidered to be described by the ACDM model, with zero spatial curvature. We assume 
Gaussian initial conditions, and therefore a vanishing bispectrum at initial time. 

During the matter dominated era and within the Hubble radius, the matter velocity 
dispersion is small. It is therefore possible to use the Newtonian limit of the Einstein's 
field equations for the perturbation of the energy tensor. By starting TRG computations 
at a redshift Zi n i = 35, and dealing only with non- linear equations for wavenumbers above 
2.10 -3 h/Mpc, one stays in the range where these approximations are valid. 
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In the TRG approach and within the ACDM paradigm, one has to renormalize both 
baryon and CDM. This could lead to a tedious system with mixed variables, since a priori 
baryons have non-negligible pressure (sound speed). However, pressure effects appears at 
very small scales (in any case, above 10 /i/Mpc), and it is safe to approximate baryons 
as pressureless and collisionless. Then, both baryon and CDM follow the same non-linear 
equations, and share the same initial condition provided that the starting redshift is chosen 
sufficiently far from decoupling (which is the case with Zi n i = 35). They can be described 
as a single perturbation field, hereafter denoted as the matter field with subscript m, and 
all quantities X m are defined as Q m X m = flbX^ + fl c X c . 

The non-linear collisionless Boltzmann equation of a self-gravitating fluid is the mas- 
ter equation describing structure formation, and being solved in N-body simulations. In 
Appendix A, we recall the basic steps and assumptions allowing to reduce the Boltzmann 
equation to a coupled system of non-linear equations in Fourier space: namely, the continu- 
ity and Euler equation for density perturbations 5 m (k, r) and velocity divergence 6> m (k, r). 

In order to solve these equations, the most common approach consists in expanding 
5 m and 9 m in perturbations, and solve the system at a given order. The TRG relies on a 
different expansion, in connected n-point correlation functions. The two approaches start 
however from the same equations. One needs to introduce simplified notations similar to 
those in j|, and to define a doublet ip a (a=l,2): 

(M) \ 

{wMJ \-e m (k, v )/n) ' 1 • ; 

where r] = log(a/ai n ;). The role of the exponential term consists in factorizing out most 
of the time evolution. Indeed, the variables <pi(k., rj) are constant during matter domina- 
tion and in the linear regime; their subsequent slow variation is easier to catch numeri- 
cally than the fast variation of (S m ,9 m ). The non-linear continuity and Euler equations 
(eqs. ( A.16| , A~T7| ) of Appendix A) can then be cast in a compact form: 



d v ip a (k, rf) = -Q a6 (k, r/)y> 6 (k, V) + e ?? 7af> c (k, -p, -q)^(p, f?)<£c(q, v), (3-2) 

where we follow the Einstein convention for repeated indices (sum) and momentum (integral 
on p and q). The vertex functions 7 a fe c (k, p,q) (a,b,c=l,2), encoding the non-linearity, are 
defined through: 

7m(k,p,q) = ^<5z)(k + p + q)a(p,q), (3.3) 

7222(k,p,q) =<5 D (k + p + q)/3(p,q), (3.4) 
7n 2 (k,p,q) = 7i2i(k,q,p), (3.5) 



and all other components vanish. The definition of the kernels a and /3 is given in eqs. ( A.l g| ) 
of Appendix [A]. The ft matrix, accounting for the linear evolution, reads 



1 -1 



(3.6) 
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including a factor B(k, rj) = Yli^m ?X^fc (* n ACDM, i runs only on photon and neutrino 
species) which appeared thought the Poisson equation. Although we are interested in 
following the matter perturbations non-linearly, we do not want to abandon altogether 
other species which might play an important role at the percent precision. The equations 
presented in this section include this effect, while those solved by the code do not. This issue 
will be discussed in more details in sections || and [7|. A nice property of this formulation is 
that everything one would like to test about new cosmologies is entirely encoded in the 
matrix, while the non-linear vertices are model- independent. 

The equations of evolution for correlation functions of ip a can be derived by applying 
several times equation (p.2[), omitting arguments for clarity: 

d v < <p a <Pb > =~ Mac < <fc<Pb > -^bc < fafc > 

+ e V Jacd < PcPdPb > +e V 7bcd < VaVcfd >, (3.7) 
9 V < PatPbVc > = - £lad < PdfbVc > ~^bd < VaVdVc > ~^cd < ^afbfd > 

+ e r? 7a (ie < ip d ip e ip b ip c > +e v j bde < ip a ip d ip e ip c > 

+ e r? 7 C( fe < (f a tf b (f d fe >, (3.8) 

Or, < PatPbtfcfd > =■■■ 



The power spectrum P ab , bispectrum B abc and trispectrum Q abc d of the doublet components 
are defined as: 



< </> a (k,r/)</? 6 (q,r/) > = 5 D (k + q)P ab (k,ri) 

< Va(k, r/)</? 6 (q, 7/)</? c (p, T})> = 5 D (k + q + p)B abc (k, q, p; rj) 

< ¥?a(k, r?)^ 6 (q, rj)(p c (p, rj)ip d (j, rj) > = 

[8 D (k + q)MP + j)Pab{k, r])P cd {p, rj) 

+ Mk + p)Mq + j)P ac (k,7/)P M (q, V) 
+ S D (k + j)<J D (q + p)P ad {k, 7?)P 6c (q, rj) 

+8 D (k + p + q + j)Qabcd(k, q, p, j; rj)] (3.9) 



One has to note that P ab will differ from the usual matter /velocity power spectrum by a 
factor e 2ri and some factors of 7-L(t) in the case of P$q and Pgg. To close this system of 
equations, the approximation proposed by M. Pietroni consists in setting Q a bcd to zero. 
Hence, the TRG does not consist in truncating at a given order in perturbations or loops, 
but at a given order in connected correlation functions. Assuming Q a bcd = should be valid 
up to some time in the non-linear evolution. This truncation still allows for a departure 
from gaussianity (the bispectrum is not set to zero at any time), and for evolving models 

one 



with a non-zero primordial bispectrum. By putting together (|3.9| ) and (3/7, 3 



finds a closed system of equations, corresponding to the TRG equations truncated at the 
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trispectrum order: 



cLP a6 (k, 77) = - fi ac (k, i])P cb (k, 77) - ft& c (k, r])P ac (k, rf) 



e v J d 3 q [7acd(k, -q, q - k)B bcd (k, -q, q - k; rf) 

+B acd (k, -q, q - k; r])j bcd (k, -q, q - k, 77)] , (3.10) 
dr,B abc (k, -q, q - k; 77) = - £l ad (k, r])B dbc (k, -q, q - k; rf) 

- ttbd{-<i, v) B adc(K -q> q - k; ??) 

- ^cd(q - k, r))B abd (k, -q, q - k; 77) 

+ 2e v [7ad e (k, -q, q - k)P db (q, r])P ec (k - q, rf) 
7bde(-q, q - k, k)P dc (k - q, ??)P ea (k, 77) 
7cde(q - k, k, -q)P da (k, 7?)P eb (q, 77)] . (3.11) 

These definitions are in agreement with the Fourier transform conventions adopted in this 
paper, namely, A(k, r) = (2ir)~ 3 J d 3 xexp(— ikx)^4(x, r). The other common definition, 
actually used in CLASS, reads A(k,r) = f d 3 xexp(— ikx)A(x, r). In this case, in the RHS 



of eq. ( 3.11 ), the last term picks up one extra (2ir) factor. 



4. Numerical implementation 

4.1 Integrated form of the equations 

Our module does not rely on the full TRG equations, but on their integrated form already 
introduced in the original TRG paper ||, that we briefly summarize below. 

As a consequence of the statistical homogeneity and isotropy of the universe, the power 
spectrum P only depends on k = |k|, while the 7 and B functions that appear in ( 3.1Q ) 



under the summation sign only depend of |q|,|k — q| and the angle between k and q. Let 
us introduce a spherical coordinate system (r, 8, ip) such that k = (k, 0, 0) and q = (q, 0, ip). 
When integrating over d 3 q, the integral over ip yields a 2tr factor. Then by performing the 
variable change between 6 and p defined as p 2 = \k + q\ 2 , one can transform the measure 
to: 



277- r°° ri+ k 
d 3 q=— / qdq / pdp. (4.1) 

J\q-k\ 



k 



Finally, by analyzing the shape of the integration domain, it is possible to further simplify 
the integrand in ( |3.10| ) and get 

d v P ab (k, rf)=- Uac(k, rf)P cb (k, rf) - tt bc (k, rf)P ac (k, rf) 

4vr f°° f q 
+ eT] -r / idq / pdp[j acd (k,q,p)B bcd (k,q,p;r]) 

k Jk/2 J\q-k\ 

+B ac d(k,q,p;vhbcd(k,q,p)] , (4.2) 

where "y a bc(k,q,p) = 7 a & c (k, q, p)| p =_(k+q) • m a numerical implementation, the bispec- 
trum, being a function of three variables and three indices, is a large object. To reduce 
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the dimensionality of the problem, it was already noticed by M. Pietroni that as long as 
the scale-dependence of the O matrix can be neglected, we can obtain a closed system in 
which the integrated quantities replacing the bispectrum read 

poo rq ^ 

Iacd,bef(k) = qdq pdp- [>y acd (k,q,p)B e f g (k,q,p) + (q «-> p)] . (4.3) 

Jk/2 J\q-k\ 1 

According to eq. ( |3.11 ), the time evolution of I is given by 
dr]Iacd,bef{k) = ^bglacd,gef{k) ^eg^acd,bgf(k) ^fg^acd,beg 

(k)+2e"A acdibef (k), (4.4) 

with the ^4's being given by 

A ac d,bef( k )= dqql dp p-{'j acd (k,q,p)[<y bgh (k,q,p)Pg e (q)P h f(p)+ 

Jk/2 J\q-k\ 1 

7egh{q,P,k)P gf (p)P hb (k) +~ffgh(p,k,q)P gb (k)P he (q)} + . (4.5) 

In eq. (|4.2| ), all integrals have been absorbed in the definition of /: 

d v P a b{k) = -ft ac P cb (k) - VL bc P ac (k) + e v — [I acd , bcd (k) + I bcd ,acd{k)). (4.6) 

Finally, we can perform a rotation from the (q,p) basis into a new (x,y) basis in order to 
obtain separable integration ranges: 



oo rq 

dq q dp p [F(k, q,p) + (q -H- p)] 

k/2 J\q-k\ 



oo rk/V2 x 2 _ y 2 

dx / dy 



_,, x + y x — y . , . 



(4.7) 



'fe/V2 

where F is a generic notation for the integrands appearing in eq. ( [4.51) . 

We implemented this integrated approach in our module as a first step (which is al- 
ready far from trivial from the numerical point of view). Having checked this level of 
approximation, we aim to go back to the full set of equations in the future. This would 
allow to compute the full bispectrum evolution, starting from Gaussian or non-Gaussian 
initial conditions. It would also allow to treat consistently the case of a matrix 0, with 
a significant scale-dependence, and to perform a more robust test of the approximation 
used in [|l^, pl[j , in which the above method was employed in the presence of massive 
neutrinos or dark energy perturbations. 



By analyzing the symmetries of eq. (|4.5|), (4.3) and (fO|), (|374|), one can reduce the 
number of independent I ac d,bef(k) functions from 64 to 14: namely, the set (121, def) 
with (def) = (111), (211), (121), (112), (122), (212), (221), (222) and the set (222, ghi) with 
(ghi) = (111), (211), (121), (122), (212), (222). Hence, we must solve a system of 17 differen- 
tial equations (3 for P ab and 14 for I ac d,bef)- The only time-consuming step is the evaluation 
of the 14 two-dimensional integrals of eq. ([T^) at each time step. 

4.2 The three modes 

Our numerical module can be tuned in such a way to solve either the above equations, or 
some limit of these equations, in order to reproduce different approaches. 



-9- 



4.2.1 Linear mode 



First, by setting all the A coefficients to zero, we can integrate the linear part of equations 



fl4.4| , |4,q) and try to reproduce the linear evolution of the power spectrum between z\ n i 
and today. Comparing these results with very accurate ones produced directly by the 
Boltzmann code is actually a useful way to estimate the precision of the time integration 



algorithm inside the TRG module, as we shall see in section 4.3. This calculation is 



performed when the CLASS input parameter "non linear" is set to "test-linear". 



4.2.2 TRG mode 

The full solution of the TRG equations, computed when the CLASS input parameter "non 
linear" is set to "trg", is of course much more time-consuming, since 14 independent A 
terms must be computed at each time step. The code has been parallelized with OpenMP. 
Its running time is of the order of 70 minutes on a single 2.3GHz core for the set of default 
accuracy parameters described in the next subsections and in Appendix [B], and on a multi- 
core machine the running time scales almost linearly (i.e. 5 minutes on 14 cores, which is 
an optimal choice due to the structure of the equations). 

4.2.3 Dynamical 1-Loop (D1L) mode 

Next, by computing the integrals of eq. ([0|) using the linear power spectra at each time, 
one expects to reproduce standard results from one-loop perturbation theory. In this 
paper, we will call this method D1L for Dynamical 1-Loop. As explained in ||, in an 
Einstein-De-Sitter universe, the D1L equations can be integrated over time analytically, 
and then reduced to eq. ( |2.1[ ). In realistic cosmological models, the linear spectra P a b are 
time dependent, and the 14 A terms need to be evaluated at each time, leading to an 
algorithm essentially as time-consuming as the full TRG one. However, as long as we stick 
to ACDM models, we can use the fact that the linear spectra evolution is accounted by 
scale-independent growth factors, in such a way that P\ 1 i eax (k,r]) oc g 2 (i]), P^f 681 °£ f 2 , 
and P^ 2 near oc gf. So, the time dependence of the integrals of eq. (fD^) computed with the 
linear spectra factors out, leading to a very fast algorithm. 

In practise, before starting the computation, we get g(rj) and f{rj) at each time-step 
using other CLASS modules. We compute the integrals of eq. ( |4.5| ) at initial time 7]^. Then, 
at each new time step, we source the TRG equations with the same integrals rescaled by 
the appropriate number of factors (g/gim) and (f/fim). This calculation is performed in 
the publicly released module when the CLASS input parameter "non linear" is set to "one- 
loop". It provides one-loop results for ACDM models with no approximation regarding the 
time evolution of perturbations, while being nearly as fast as the S1L method. We will see 
in section [6| that it greatly improves over S1L predictions. 

4.3 Time integration 

The time integration is performed with a predictor-corrector method. Starting from a time 
rj at which all quantities are known, one computes the new P and I functions at the step 
rj + Srj/2 with the standard Euler method, as well as their derivatives in this point. The 
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values of P and I at rj + Srj are then derived with the Euler method, using however the 
derivatives at the mid-point. With respect to a basic Euler method, this allows to take 
many less steps for the same level of accuracy (typically, we use between 50 and 100 time 
steps to evolve from z = 35 to the present time). 

With a predictor-corrector method on 100 time steps, and at all redshifts and wave- 
numbers of interest, we find that the linear spectrum calculated by the TRG module agrees 
at the 0.05% level with the one calculated by CLASS (with many more time steps and a 
sophisticated integrator). With 50 time steps, the error is still as low as 0.08%, allowing us 
to keep this choice as a default setting in the released version of the code. The convergence 
has also been controlled for the one-loop computation: with only 50 steps and a predictor- 
corrector method, one reaches the same 0.1% level accuracy as with lO'OOO steps and a 
standard Euler method. 

Convergence of the time integrator - linear mode 

0.0009 
0.0008 
0.0007 
\ 0.0006 

ro 

I 0.0005 

1? 0.0004 

I 0.0003 

Ql 

0.0002 
0.0001 


0.5 1 1.5 2 2.5 3 3.5 

r|=log(a/a ini ) 

Figure 2: Convergence proof of the linear time-evolution in the TRG module w.r.t. the linear 
evolution in CLASS. 

4.4 Momentum integration 

Since in the TRG module, most of the time is spent in the computation of A functions, an 
obvious way to speed up the code would be to make a cut in the domain of integration. 
This domain, defined in eq. fl4.7| ), has a rectangular shape (x > k/y/2, < y < k/^/2). 

After noticing that the main contribution to the integral comes from the region near 
the x = k/y/2 and y = k/y/2 boundaries, and that the integrand exhibits an approximate 
symmetry around the axis y = V2k — x, we first tried to integrate on L-shaped bands 
starting from the top left corner in (x, y) space. One could hope that after integrating 
over a few such L's, the integration would converge and could be stopped. This does not 
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work, and we noticed that the tails of the integrand over the whole domain contribute 
enough to create departures from the exact result when implementing a cut. Since no good 
trade-off was found, we kept the whole domain. We use however logarithmic integration 
steps in the (x,y) space, defined in such way to (i) increase the resolution near the top 
left corner in (x,y) space and (ii) preserve the approximate axial symmetry. Moreover, the 
value of these logarithmic step sizes is modulated as a function of the wavenumber k. The 
key parameter controlling the trade between precision and overall speed of the module is 
called logstepx_min. This parameter provides a lower bound on the logarithmic step used 
in both x and y directions, for any wavenumber k. We show in figure || that setting this 
parameter to 1.06 is sufficient for the final result to converge at the 0.6% level in the one- 
loop case (at z = 1 and for k < 0.5/i/Mpc), and 0.3% level in the TRG case. Performing 
the test at z = gives similar results. 



one-loop mode, z=1 



trg mode, z=1 




0.2 0.3 

k (h/Mpc) 




0.2 0.3 

k (h/Mpc) 



Figure 3: Convergence proof for the integration in momentum space in the the one-loop (left) 
and TRG (right) modes. The reference model has logstepx_min=1.001 in the one-loop case, or 
logstepx_min=1.02 in the TRG case. 

In all integrals, we define some ultra-violet and infra-red cuts-off (fc m in 5 &max) besides 
which the power spectra are assumed to vanish. Moreover, below a scale ki, we assume 
that the evolution of the spectrum remains strictly linear at all times: this means that for 
k < ki, we set the A coefficients to zero, in order to speed up the calculation. We checked 
that the results are well converged with the choices k m \ a = 10 _4 /i/Mpc, k^ = 10 _3 /i/Mpc, 
U = 10 3 VMpc. 

4.5 Numerical instabilities 

The main challenge in any implementation of the TRG algorithm is to keep numerical noise 
under control in the ultra-violet limit. Indeed, a small error in the tail of the power spectrum 
would be amplified by the recursion over time, and lead to exponential divergences. By 
the same phenomenon, a slightly non-smooth initial spectrum would produce non physical 
divergences. The smoothness of the spectrum is very difficult to control since it is passed 
in the form of tabulated values. 
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We checked that the divergences and wiggles observed in our first, straightforward 
implementation of the algorithm were only numerical artifacts. Indeed, we varied the size 
of time and wavenumber steps, and observed that these features propagated at a rate 
defined by the number of discrete steps, rather than by physical time. 

In order to eliminate these numerical artifacts, we tried to implement several smooth- 
ing algorithms and several exponential cut-off functions in the small-scale linear power 
spectrum. Every time, some feature would finally propagate. By the very nature of the 
non-linear equations, as time goes by, any small error is amplified in huge proportions. 

We solved this issue with an approach which is physically justified by the fact that 
during structure formation, mode coupling leads to a transfer of power from large scales 
to small scales, and not in the other direction. At each new time step, numerical artefact's 
unavoidably affect the last four points in the list of discrete k values (two points in the 
computation of the time derivative, plus two points in the calculation of the A integrals). 
Hence, one can systematically throw away these last four points. In practice, this method 
(that we call "double escape") consists in integrating the TRG equation only in a range 
kh < k < &de(?7)- The wavenumber /cde(??) (where DE stands for double escape) coincides 
with k max at initial time, and then decreases by discarding four more points in A;-space at 
each half time step Srj/2. Since a given wavenumber is impacted by the evolution of smaller 
wavenumbers only, our method cannot affect the final results. We checked it explicitly by 
changing the step sizes in k space. Such a reduction means that we decrease /cde(^) 
differently as a function of time, but the final results remain totally invariant. Simply, 
smaller step sizes allow to get results at a larger /cde at z = 0, starting from the same 
initial k ma , x . In the released version, the default time steps are such that /cde ~ 0.25/i/Mpc 
today. When the user asks for non-linear spectra at a given redshift, they are only written 
in output files up to the corresponding value of &de(??)- 

Note that when computing the A integrals, the power spectra still need to be evaluated 
between &e>e(??) and /c max , but the final results have a very weak dependence on P a b{k) in 
this range. We use a log-linear extrapolation of P a b{k) besides &de(?7), based on an estimate 
of the logarithmic derivative in k = /cde(?7)- We actually tried alternative extrapolation 
schemes (e.g. with constant second derivative in log-log space) and found that the final 
result is independent on the scheme (again, as a consequence of the fact that a wavenumber 
is only affected by its smaller neighbors) . 

We conclude that this "double escape" is a simple and robust way to keep numerical 
instabilities under control, leading to a "loss of information" on the smallest wavelengths, 
but introducing no bias in the final results. 

Note that these problems arise when dealing with the full TRG equations. In its one- 
loop limit D1L, which will be seen below to provide excellent results, the integration over 
time is perfectly stable, and there is no need to throw points away. 

5. Self-consistency checks in the exact EdS limit 

The main purpose of this section is to test the agreement between the S1L and D1L 
methods. This check is more straightforward in a model without non-trivial growth factors 
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Figure 4: In an ideal EdS limit and at four different redshifts, comparison between the non-linear 
spectra computed with the standard one-loop S1L, dynamical 1-loop D1L and TRG methods. For 
the D1L and TRG case, we start the computation from exactly linear initial conditions at z; n ; = 35 
(with vanishing bispectrum), causing a residual difference with respect to the S1L method. 

related to Vl m begin smaller than one, i.e. in the Einstein-De Sitter limit. Hence, for the 
results of this section, we modified our TRG/D1L module in order to enforce a perfect EdS 
evolution. This amounts in: 

• replacing the actual matrix elements by (3/2, —3/2), 

• for D1L, replacing the exact growth functions by g(rj) = 1, f(rj) = 1. 

The issue of initial conditions is crucial. The simplest choice consists in starting from 
exactly linear perturbations at z = z- m i, i.e. from a Pu(k) given by the linear density 
power spectrum, from P22(k) = Pu(k) = Pu(k), and from I ac d,bef{k) = 0. Starting from 
such conditions at z- m \ = 35, we obtain the 1-loop density power spectra named D1L in 
figure ||[ They are actually smaller than S1L results, at all redshift and by several percents 
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Figure 5: Same as previous figure, although for the D1L and TRG case we now start from an 
initial condition / = 2A, which amounts in sharing the same assumptions as the S1L calculation. 

(at k = 0.25 /iMpc -1 they predict non-linear spectra differing respectively by 1.2, 1.5, 2, 5% 
at z = 3,2,1,0). This discrepancy actually decreases when we start the integration at 
a higher redshift. This suggests that the difference arises from neglecting any non-linear 
evolution before Z[ n {. Indeed, the S1L method assumes that from z — > oo till today, 
the modes have always evolved according to the newtonian equations in the EdS limit, 
with a linear growth factor proportional to a. This is of course an approximation, but in 
order to do a consistent comparison of the methods, we should also integrate the D1L and 
TRG equations starting from initial conditions reflecting the 1-loop evolution in the range 
z E [oo,Zi n i], instead of linear initial conditions. 

This amounts in taking the initial P a b{k) to be the one-loop spectra at rj = 0, and in 
computing I ac d,bef(k) at rj = at leading order in perturbation theory. However, one-loop 
corrections to the linear P a b(k) at z ~ 35 appear to be completely negligible in the range of 
interest (at least until k ~ 1 ftMpc -1 ), so sticking to the linear spectra makes no difference 
in practice (we checked this explicitely) . This is not true as far as the initial bispectra 
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are concerned: instead of B = 0, we should use the tree-level bispectrum expression ||, 
computed at z- m \ out of the linear power spectra. When plugged into the expression of 
Iacd,bef(k), the tree-level bispectrum leads to a series of terms that we computed explicitely. 
We checked numerically that all terms are negligible with respect to the leading one, namely 
/ = 2 A: so, for simplicity, we will approximate the tree- level expression of / as 2 A. A more 
intuitive derivation of this result follows from observing equation (^4|). At very early time, 
the bispectrum is driven away from zero by the source Axed,&e/) so that 

d v Iacd,bef( k ) - 2eVA acd,bef( k ) ■ 

Since at leading order in perturbation theory A ac d,bef{k) is constant in time, this approxi- 
mate evolution equation is solved by I a cd,bef{k) = 2e v A acc i t bef (k) , which corresponds to the 
initial condition I ac d,bef( k ) = 2A acdMf (k) at r] = 0. 

When adopting such initial conditions, the D1L results nicely coincide with S1L predic- 
tions, as can be seen in figure [|. This was of course expected on a purely analytical basis, 
but the agreement proves that the two methods are correctly implemented numerically. 

The most important result of this paper is already visible in Figures || and ||[ on 
the scales displayed in the figures, the TRG curves are always very close to their D1L 
counterpart with the same initial conditions, at least for z > 2. At z = 1, differences 
at the level of half-a-percent start to appear. At z = one can see clear differences of 
1 or 2% on the scales corresponding to BAO dips: the TRG smoothes BAOs more than 
1-loop approaches. However, this difference between D1L and TRG results remains small, 
especially in regard of the large increase in computing time in the TRG case. Since in a 
realistic scenario the growth of structure is suppressed with respect to the EdS case, this 
conclusion is likely to hold even better in ACDM, as we shall see below. 



6. Results for ACDM and comparison with N-body simulations 

We finally ran our code for a ACDM model, in order to check that the findings of the 
previous section remain true, namely: (i) the fact that the D1L method with / = 2 A 
matches S1L predictions (despite the fact that D1L treats in a more exact way the impact 
of Q m < 1 at early time due to the small radiation density, and the consequence of late A 
domination for the time evolution of the density and velocity linear growth factors), and 
(ii) the fact that TRG results do not improve significantly over one- loop calculations. 

Also, in order to evaluate the precision of each method, we wish to compare semi- 
analytic predictions with very accurate N-body simulations, leading to a density power 
spectrum with a good resolution in /c-space, and a small sampling variance on mildly non- 
linear scales, hopefully quantified by an error bar. Such power spectra can only be obtained 
by running several simulations in a very large box, and computing the mean and variance 
of the resulting collection of spectra. Simulations with these characteristics have been 
presented in several papers including Refs. [||, ||, 16]. In Refs. [||, 16|, the initial linear 



spectrum was inferred from analytical fitting formulas, which are not accurate at the per- 
cent level. We prefer to use the recently published simulations of Sato et al. ||, based 
on an initial spectrum computed with CAMB. These authors performed a large number 
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Figure 6: Matter power spectrum divided by a smooth one [HJ for z = 3,2,1,0.5 obtained with 
five different methods: N-body simulations TRG method, standard one-loop S1L, dynamical 
one-loop D1L, and finally HALDFIT [Q. Here we employed non-physical initial conditions / = 2 A in 
order to check the agreement between D1L and S1L results, but we know that these results cannot 
be trusted. 



of high-quality simulations, leading to sub-percent level statistical errors, and to a very 
good matching between linear and N-body spectra for wavenumbers in the range [0.05 — 
0.1]/i/Mpc. 

We show in figures ^, a comparison between those results and the non-linear spectra 
obtained with the one-loop (S1L and D1L), TRG and HAL0FIT semi-analytical methods. 
All spectra are computed for a model with h = 0.71, SI a = 0.735, fi& = 0.0448, iV fj = 3.04, 
T cmb = 2.726K, n s = 0.963 and a s = 0.80. 

Like in the previous section, we have to face the crucial issue of initial conditions. 
We know that starting from the tree-level bispectrum, equivalent to / = 2A, amounts in 
doing the same assumption as in the S1L calculation: it assumes that the bispectrum has 
grown like in an ever-matter dominated universe, containing a single presureless species 
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Figure 7: These figures contain our main results. The matter power spectrum divided by a 
smooth one for z = 3,2,1,0.5 obtained with five different methods: N-body simulations (§), 
TRG method, standard one-loop S1L, dynamical one-loop D1L, and finally HALOFIT 0|. In contrast 
with the previous figure, we now assume linear initial conditions ar Zj n j. This assumption appears 
to provide excellent results on all displayed scales till z = 2, and still better predictions than halofit 
at z = 1,0. 



described by Newtonian equations. Actually, most recent N-body simulations (including 
the one quoted here) also share this assumption. Indeed, their initial conditions are set 
by a 2LPT (2nd-order Lagrange Perturbation Theory) algorithm, designed in such way to 
minimize transient effects. This result is achieved precisely by assuming a tree- level initial 
bispectrum at initial time. Hence, the choice I = 2A is the proper one for comparing D1L 
or TRG results both with S1L predictions and with N-body results. 

We show in figure || the D1L and TRG results obtained from such initial conditions. 
The S1L and D1L method are again in very good agreement. We could expect D1L and 
S1L to depart from each other at z < 1, when A domination leads to non-trivial growth 
factor which render the expansion of eq. (2.4) only approximate. This difference can be 
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observed at z = 0.5 on the scale k = 0.125 /i/Mpc corresponding to a maximum in the 
BAOs. However, it is only a sub-percent effect, which confirms the results of Ref. 0, 
figure 7, in which a similar test was presented. As in the previous section, the TRG results 
are very close to the S1L ones , even at z = 0.5. All three methods tend to overestimate 
non-linear corrections to the density power spectrum: this trend was already well-known 
for calculations limited to one-loop, and our result shows that TRG results make roughly 
the same error. Percent precision is nevertheless achieved until k ~ 0.2/iMpc^ 1 for z > 2, 
or until k ~ 0.14 /iMpc^ 1 at z = 1. 

For comparison, the results obtained when starting from 1 = 0, which are presented 
in figure f|], seem to be in much better agreement with N-body simulations. Note that 
earlier comparisons shown in Ref. g ||, @, (jj, 0] 

were based on the same assumptions. 
However, this comparison is not fully self-consistent, since it relies on different choices of 
initial conditions in the different methods. Hence, the better agreement observed in this 
plot is likely to be purely coincidental. 

7. Discussion 

In this work, we compared a few semi-analytical methods to the recent accurate N-body 
simulations of Sato et al. || and to standard 1-loop perturbation theory. We released these 
methods in the form of a C module, named trg . c and integrated in the Boltzmann code 
CLASS, version 1.2. The first method is the Time Renormalization Group (TRG) proposed 
by M. Pietroni. The second method, named here D1L (Dynamical 1-Loop) and based 
on the same approach, leads to nearly identical results on scales of interest, with a much 
smaller computing time. We have explained in this paper our approach for dealing with 
numerical instabilities, and for optimizing the algorithm. We have also presented some 
convergence tests proving its reliability. This release allows any cosmologist to compute 
easily an approximate non-linear spectrum (from first principles, rather than with fitting 
formulas), while until now this task was only accessible to a few specialists. 

Differences between TRG and D1L are visible at low redshift, but they remain very 
small on scales of interest. This shows that the partial resumation of diagrams beyond 
one loop in the TRG method improves one- loop results by a negligible amount, for a much 
larger computing time. Hence, the one- loop limit of the TRG equations (namely, the D1L 
scheme) is preferable in practice. 

The D1L algorithm is a fast and practical tool for obtaining one- loop results for the 
density/velocity power spectra (and tree- level results for the bispectra), even for non-trivial 
cosmological models or in presence of an arbitrary initial bispectrum. There are several 
ways to incorporate such assumptions in one-loop calculations, but in the D1L equations, 
this flexibility is present from the beginning. In this paper, we used this opportunity for 
showing the importance of assumptions concerning the initial bispectrum. In our universe, 
the non-linearity of the gravitational equations is sufficient to induce a tiny bispectrum even 
at very high redshift on cosmological scales of interest in this paper. At any given initial 
redshift, this bispectrum is small enough to be well approximated by a tree-level calculation, 
but nevertheless large enough to impact results at small redshift by several percents. This 
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observation is consistent with previous studies of initial conditions for N-body simulations. 
In particular, the 2LPT method has been developed in order to deal with transient effects 
in simulations [jT^ ], i.e. in order to remove spurious decaying modes by implementing 
initial conditions infered from second-order Lagrangian perturbation theory, including an 
initial tree-level bispectrum. However, this fact had been overlooked in the context of TRG 
calculations. We have shown that previous claims that TRG results improve over one-loop 
predictions were the consequence of neglecting this initial bispectrum. When using the 
same initial bispectrum in all methods, the difference between S1L, D1L and TRG results 
almost disappears. 

It is beyond the scope of this paper to discuss to which extent the initial conditions 
assumed in the S1L method and in N-body simulations offer a sufficiently realistic de- 
scription of the actual universe. The fact that at high redshift baryons, cold dark matter 
and radiation perturbations do not share the same transfer functions and do not obey ex- 
actly to the growing mode solution means that in principle, one should implement some 
amount of "physical transient effects" in the initial conditions of both N-body simulations 
and renormalisation algorithms. The authors of have already shown the importance 
of treating baryons separately from cdm when performing non-linear calculations. To our 
knowledge, the impact of early baryon and radiation perturbations on the calculation of 
the initial bispectrum (to be passed to N-body initial condition generators) has not been 
quantified accurately. We speculate that the D1L algorithm described in this paper could 
be a convenient tool for computing realistic high-redshift bispectra, and we leave this for 
future studies. 
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A. Prom Boltzmann to continuity and Euler equation 

Let's start by considering a large set of particles only interacting through gravitation. In 
order to lighten slightly this presentation, let us omit the m subscript for the moment. For 
each particle (CDM or baryon) the equation of motion is given in terms of its proper time 
t, velocity v and position r: 

where (ft is the Newtonian potential induced by the local mass density p(r), 

<ft{r) = G f dV-^-^. (A.2) 
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This set of particles is however sitting in an expanding universe. The gravitational collapse 
consists in studying departures from the homogeneous Hubble expansion. It is then useful 
to redefine the variables of position and momentum with respect to the comoving ones. 
The physical space coordinates r relate to the comoving space coordinates x through r = 
a(r)x. We can define the density perturbation 5(x,t), the peculiar velocity u(x, r) and 
the cosmological gravitational potential $(x, r) by 

P (x,t) = p(t)[1 + <5(x,t)], (A.3) 
v(x,t) = ^x + u(x,t), (A.4) 

0(x,r)^-^x 2 + cl>(x,r). (A.5) 

The Einstein equations lead in the sub-Hubble, weak gravitational field limit to the Poisson 
equation: 

V 2 $(x,r) = -|?£ 2 ^n < (r)<5 i (x,r). (A.6) 



Defining the momentum p = amu, eq. ( |A.1| ) now reads: 

= -arnV$(x). (A.7) 

dr 

One finally write the collisionless Boltzmann equation for the phase-space density /(x, p, r) 
in its full form: 

df df p df 

— = — H V/-amV$-— . (A. 8) 

dr or ma ap 

The non-linearity of this equation is induced by gravitational back-reaction, described by 
the Poisson equation. This equation is obviously hard to solve due to its high dimension- 
ality, which can however be reduced by taking the momenta of the distribution /(x, p, t), 
and by truncating such an expansion at some order, if this can be physically justified. Let 
us define the first three momenta as: 

yd 3 p/(x,p,r) =p(x,r), (A.9) 

d 3 p -?-/(x, p, r) ee p(x, r)u(x, r), (A.10) 

am 

/ d 3 p -^-4/(x,P,t) EE p(x,T)UiUj(x,T) + (Tij(x,T). (A. 11) 

7 a z m z 

Taking the zeroth moment of eq. ( |A.8j ) leads to the continuity equation: 

^2ll + V • [1 + 5(x, r)]u(x, r) = 0, (A.12) 
while the first moment gives (after subtracting u(x, r) times ( A.12| )) the Euler equation: 
dUl( ^ T) + H(r) Ui (x, r) + (u(x, r) • Vu(x, r)) l = -V^(x, r) - iv^-. (A.13) 
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One usually takes the anisotropic pressure cr^- to be equal to zero. This approximation is 
excellent not just for baryons (they interact, so they are considered as a perfect fluid with 
no shear), but also for cold dark matter. Indeed, CDM particles have no interactions but 
a tiny velocity dispersion, leading to negligible anisotropic pressure on scales where CDM 
can be represented by a single-flow distribution. In the range of linear and quasi-linear 
scales considered here, this approximation is valid. 

With such considerations, we can safely use from now on the S m notation for both 
baryons and CDM. The Poisson equation (|A.6| ) can be cast into a more convenient form: 

V 2 $(x,t) = -ln 2 n m {r)5 m {^r) I 1+ £ , (A.14) 

where in i stands for photons and neutrinos (relativistic or not). It can be seen from 



eq. (|A.13|) that the curl part of the velocity field decouples and is suppressed by the universe 
expansion. The equations can then be reformulated in terms of the density perturbation 
and the velocity divergence 9 = V • u. Furthermore, to have a better understanding of the 
mode coupling induced by the non-linearity, we transform the equation into Fourier space, 
with the following convention: 

A(k, r ) = * | d 3 xe-^(x, r). (A.15) 

In Fourier space, the equations governing the matter field reduce to 
<9<5 m (k,r) 



dr 

d9 m {k, r 



+ 9 m (k,T) = - J d 3 qd 3 pfe(k-p-q)a(p,q)^ m (p,r)A m (q,r), (A.16) 



Or 



+ n(T)9 m (k, T ) 



+ ^( r ) 2 E^( T )^( k ' T ) = "/ d 3 qd 3 pfe(k-p-q)/3(p,q)a m (p,r)0 m (q,r), (A.17) 
where the kernels a and (3 encoding the mode coupling are defined as: 



(p + q) • p (p + q) 2 (p • q) , 

a(p,q) = 2 — ' /3p,q) = ttt^ • A - 18 



B. Structure of the trg.c module in CLASS 

As explained in ||, CLASS is organized in eleven modules. The role of the ninth module 
nonlinear. c is to evaluate the non-linear power spectra according to the method chosen 
by the user. In the current version CLASS vl.l, the nonlinear, c module can optionally 
call the trg.c module described in this paper, in one out of three modes (non linear = 



test-linear, one-loop or trg) already described in subsection L2), and with a set of 
precision parameters. After the execution of the trg . c module, the code is ready to write 
in output files the non-linear density, velocity and cross spectra at different redshifts chosen 
by the user. 
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The trg . c module consists of several functions. The most important one, trg_init () , 
is called from the nonlinear, c module. Its goal is to compute the non-linear power 
spectrum from a given starting redshift and linear power spectrum previously computed 
by the spectra. c module. 

The trg_init() routine first defines relevant step sizes in ij space and in k space. For 
the latter, it calls the two functions trg_logstepl_k() and trg_logstep2_k() . The k 
steps are defined in such way to keep a high resolution in the baryon-oscillation zone, and 
a reasonably small number of steps on other scales. Since these steps are relevant for the 
double escape procedure defined previously, one might want to be cautious while modifying 
them. 

The initial density spectrum P\\ is taken directly from the spectra, c module. The 
initial velocity and cross spectra P22 and P\% are computed by evaluating a finite differ- 
ence between P\\ at r)i n i ± e. At the starting redshift, three-points correlating functions 
(bispectra) are assumed to be zero, so all J's defined in eq. ( fl.3[ ) are initially zero. 

The integration of the TRG equations (full equations in trg mode, simplified equations 
in the other two modes) is then performed. At each step in the trg mode, or only once 
in the one-loop mode, trg_init() calls trg_integrate_xy_at_eta to perform explicitly 
the integrals in (x, y) space and find each of the 14 A factors. This last function evaluates 
each of the 14 integrands at each (x,y) point by calling the function trg_A_arg_trg() (or 
trg_A_arg_one_loop ( ) ) . 

After the execution of the trg . c module, the non-linear spectra are stored in a structure 
nonlinear associated with the nonlinear, c module. They can be evaluated at any value 
of k and z by calling the function nonlinear_pk_at_k_and_z() . Instead, all values at a 
given z are returned by the function nonlinear_pk_at_z () . The output module output . c 
actually calls this last function in order to write the final results in separate files for the 
density, velocity and cross spectra. 
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